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Within  the  context  of  discrete  formulations  arising  from 
a  finite  element  approximation,  explicit  expressions  for  the 
frictional  consistent  contact  tangent  stiffness  and  residual 
are  derived  from  variational  equations  by  using  a  consistent 
linearization  procedure  for  both  the  sliding  and  adhesion 
phases.  The  consistent  tangent  operator  is  always  non- 
symmetric  for  the  case  of  frictional  sliding  owing  to  the 
nature  of  the  Coulomb's  friction  law  employed. 

For  two-dimensional  applications,  a  three-node  contact 
element  is  employed  in  the  finite  element  discretization. 
Numerical  examples  are  also  presented  that  illustrate  the 
performance  of  the  proposed  formulation. 
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1.  Introduction 

Frictional  stick-slip  contact  phenomena  constitute  important  aspects  of  real  engineering 
applications.  The  micromechanics  of  friction  for  metallic  surfaces  in  contact  has  drawn 
significant  attention  in  the  literature  during  the  past  several  decades.  Basically,  the  contact 
surfaces  are  not  smooth  planes  but  rough,  uneven  surfaces  composed  of  many  asperities. 
Microstructurally,  the  asperities  experience  local  yielding  and  fracture.  Moreover,  they  deter¬ 
mine  the  real  area  of  contact  as  opposed  to  the  nominal  area  of  contact.  For  a  review  on  the 
physical  aspects  of  friction,  see  for  example  Tabor  [1981],  Oden  &  Martins  [1985]  and  Oden 
&  Pires  [  1 984], 

For  a  finite  element  solution  frictionless  contact  problems,  a  perturbed  Lagrangian  and 
an  augmented  Lagrangian  formulations  have  been  proposed  by  Simo.  Wriggers  &  Taylor 
[1985],  and  Landers  &  Taylor  [1985],  respectively.  On  the  other  hand,  Duvaut  &  Lions 
[1976],  Oden  &  Martins  [1985],  and  Oden  &  Lin  [1986]  have  proposed  some  regularization 
techniques  of  the  Coulomb’s  law  of  friction  for  numerical  solution  of  dynamic  frictional  con¬ 
tact  phenomena. 

In  this  study,  a  perturbed  Lagrangian- based  formulation  is  proposed  for  the  finite  ele¬ 
ment  solution  of  fully  nonlinear  frictional  contact  problems.  The  stick-slip  contact 
phenomena  is  accommodated  by  means  of  a  nonlinear  variational  formulation.  In  view  of  an 
analogy  between  the  Coulomb’s  law  of  friction  for  the  stick-slip  motion  and  the  yield  criterion 
for  classical  elastoplasticity  (see,  e.g.,  Michalowski  &  Mroz  [1978]).  a  two-step  operator  split¬ 
ting  methodology  is  employed. 

In  the  current  literature,  the  modification  of  the  tangent  stiffness  accounting  for  the  con¬ 
tribution  of  frictional  contact  often  takes  the  form  of  symmetric  and/or  non-svmmetric  rank- 
one  updates  inherited  from  the  linear  theory  (see,  e  g.,  Oden  &  Martins  [1985]  and  Flughes  et 
al.  [1976]).  In  the  finite  element  solution  of  geometrically  ronlinear  frictional  contact  prob¬ 
lems,  however,  such  simple  procedures  are  no  longer  adequate.  In  the  event  of  frictionless 
contact,  a  consistent  tangent  operator  has  been  obtained  by  Wriggers  &  Simo  [1985]. 

Within  the  context  of  finite  element  discrete  approximation,  explicit  expressions  for  the 
frictional  contact  tangent  stiffness  and  the  residual  are  derived  in  this  paper  from  variational 
equations  by  using  a  consistent  linearisation  procedure  for  both  the  sliding  and  adhesion 
phases.  It  is  shown  that  for  the  case  of  frictional  stick  the  consistent  tangent  operator  is  sym¬ 
metric  only  when  numerical  convergence  is  achieved.  On  the  other  hand,  the  consistent 
tangent  operator  is  always  non-symmetric  for  the  case  of  frictional  sliding  owing  to  the  nature 
of  the  Coulomb's  friction  law  employed.  Not  surprisingly,  these  expressions  degenerate  to  the 
classical  rank-one  corrections  of  the  stiffness  matrices  in  the  limiting  case  of  infinitesimal 
deformations.  It  is  emphasized  that,  in  the  presence  of  nonlinear  contact  kinematics,  use  of 
the  consistent  contact  tangent  stiffness  is  essential  in  preserving  the  quadratic  rate  of  asymp¬ 
totic  convergence  of  Newton’s  method. 

For  two-dimensional  applications,  a  three-node  contact  element  is  employed  in  the  finite 
element  discretization  The  perturbed  Lagrangian  -based  computational  algon  n.  is  capable 
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ot  pertorming  a  one-pass  or  two-pass  contact  slide-line  logic.  A  number  of  numerical  exam¬ 
ples  are  presented  in  Sec.  3  that  illustrate  the  performance  of  the  proposed  variational  formu¬ 
lation. 


2.  Discrete  variational  formulation 

Within  the  framework  of  the  finite  element  method,  the  governing  variational  equations 
involving  fully  nonlinear  kinematics  are  considered  in  this  section  for  both  adhesive  and  slid¬ 
ing  contact  problems.  In  what  follows,  for  simplicity,  attention  is  focused  on  the  two- 
dimensional  (planar)  geometry.  The  extension  to  three-dimensional  geometry  is  complicated 
by  geometric  considerations  only. 


2.1.  Finite  element  discretization 

Throughout  the  remaining  part  of  this  report,  we  employ  the  bilinear  isoparametric  ele¬ 
ments  for  "parent"  contacting  bodies.  Concerning  the  contact  segment  characterization,  the 
"master-slave"  slide-line  contact  logic  is  adopted  (see  Hallquist  (1983]).  In  particular,  a  three- 
node  contact  element,  consisting  of  two  "master"  nodes  and  one  "slave"  node,  is  used,  see  Fig¬ 
ure  1.  With  reference  to  Figure  1,  the  tangent  and  normal  vectors  are  defined  as  follows: 


X*> 

1*: 

-  x, 

-  X,  1 

(2.1) 

n  =  e3 

x  t  . 

(2.2) 

where  e?  denotes  the  unit  base  vector  normal  to  the  plane  of  the  three-node  element  and 
Xi  =  X|  +  U|,  x:  =  X;  +  u;  signify  the  current  positions  of  the  master  nodes  (X,  ,  X;  for  refer¬ 
ence  coordinates  and  u,  ,  u:  for  current  nodal  dispacements).  In  addition,  we  define  the 
current  "surface  coordinate"  a  as  follows 


in  which  x,  =  X,  -  u(  denotes  the  current  position  of  the  slave  node.  The  normal  and  tangen¬ 
tial  gaps  (penetrations)  associated  with  a  typical  three-node  element  are  defined  as 

=  (xs  -  x,)  •  n  (2.4) 

g,  =  (xt  -  X,)  •  t  -  a"  I  X;  -  x,  I  .  (2.5) 


where  a  '  is  the  (old)  surface  coordinate  at  the  last  time  step  (known).  The  variations  (incre¬ 
ments)  of  the  normal  and  tangent  vectors  due  to  nonlinear  kinematics  can  be  shown  to  be  (see 
Wnggers  &  Simo  [  1 985]) 


on  =  - 

1  x:  -  x,  |  U 

Q  n) 

•  (v:  -  V\) 

(2.6) 

ot  =  - 

.  in 

3  n) 

*  (*t:  -  Vi ) 

(2.7) 

V  ’V  ■  V  V*  1  V  v  ■  ■  I  *  TV  (."■ 
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where  tj  is  the  variation  (increment)  of  u  .  Furthermore,  for  convenience,  we  define  the  fol¬ 
lowing  abbreviations  (operators) 

F)  s  (  •  )s  -  ( 1  -  a )  (  •  ),  -  a(»);  (2.8a) 

FI  =  (  •  )2  -  (  •  ).  (2.8b) 

With  the  above  notations  at  hand,  we  now  give  the  variational  derivation. 


2.2.  Frictional  stick 

For  the  case  of  frictional  stick  (no-slip),  we  consider  the  following  perturbed  Lagrangian 
functional  for  bodies  in  contact: 


n„(u  ,  A,  ,  A,)  =  ri(u)  +  A^G, 


~  A-n  h.n  +  A,r  G, 

2u„ 


(2.9) 


Here  u  designates  the  vector  of  nodal  displacements,  A,  (A,)  the  vector  of  normal  (tangential) 
nodal  contact  forces,  G,  (G,)  the  vector  of  normal  (tangential)  nodal  gaps,  and  to,  (to,)  the 
normal  (tangential)  penalty  parameters.  Moreover,  Il(u)  stands  for  the  total  potential  energy 
of  the  bodies  in  contact. 

The  discrete  variational  equations  are  then  obtained  by  taking  the  variations  with 
respect  to  u  ,  A,  ,  and  A,  ,  respectively: 

5„n(u)  +  AjSuG,  +  A,r5uG,  =  0  (2.10a) 

5A,r(- —  A,  +  G,)  =  0  (2.10b) 

“n 

5A,r  (-  —  A,  +  G,)  =  0  (2.10c) 

<0, 


From  (2.10b,c),  we  obtain  that  A,  =  <o„  G,  and  A,  =  oo,  G,  as  the  normal  and  tangential 
penalty  contact  forces. 

The  variation  of  a  typical  nodal  normal  gap  gn  e  G,  takes  the  form  (see  Wriggers  & 
Simo  [1985]) 

bg„  =  VT  C,  =  [irs  -  ( 1  -  a )  JJ,  -  a  tf;]  •  n  s  ij  •  n  (2.11) 


where  c„  =  D,(i5g, )  (D  is  the  directional  derivative  operator).  Alternatively,  we  can  write 


^ n 


n  , 


(2.12) 


with  the  bar'  quantity  defined  in  (2.8a).  We  shall  give  the  matrix  representation,  within  the 
context  of  three-node  contact  elements,  of  n  later  in  this  section. 

Similarly,  the  variation  of  a  typical  nodal  tangential  gap  g,  e  G,  can  be  obtained  accord¬ 
ing  to 


bg,  =  (if.  -  I;,)  •  t  *  (x,  -  x,)  •  <5t  -  a"  h  |  x:  -  x,  | 
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=  t‘^T^r(n*’)s,,rc' 


v  s  lL.a«  =  1!  -o  -  a°)  TJ,  -  a0  7,2  , 


c,  =  D,(<5£,)  =  t°  + 


I  X2  -  X,  I 


(2.15) 


Moreover,  the  residual  vector  Rfl  and  tangent  stiffness  Kfl  associated  with  the  total  potential 
energy  of  the  contacting  bodies  simply  read 

RasD,(n(u))  (2.16) 

KssD,(Rs)  (2.17) 

In  the  case  of  inelasticity ,  Rj  and  KB  are  deduced  from  a  (Galerkin)  variational  functional  n 
involving  constitutive  relations  and  boundary  conditions. 

The  variational  equations  (2.10a,b,c)  can  now  be  stated  as 


(2.18a) 


i,  JR*  +  M  (X^c^  +  \\s)  cSs))  j  =  0 
6\U~—  A*  +  G„)  =  0 


(2.18b) 


(-  —  A,  +  G,)  =  0 


(2.18c) 


In  (2.18a),  (scalars)  X(„"  e  A„  ,  X!,"  e  A,  ,  and  M  represents  an  assembly  operation  over  all 
three-node  contact  elements  in  consideration  (5  =  total  number  of  slave  nodes  in  contact  = 
total  number  of  conditions  of  constraints).  To  apply  the  Newton's  iteration  scheme,  con¬ 
sistent  linearization  of  Eq.  (2.18a,b,c)  at  (  u  ,  A„  ,  A,  )  is  performed  and  leads  to 

Kfl  +  14  [K'„s)  +  K',s)]  14  c<f>  U  c<" 


’  nT  .  i>A‘n  ,  5A,7 


14  c!"r 


s 

i4  c'rr 

I-I 


J  -  I  s  ■  1 

I  0 

“>n 


Rfl  -  14  (  X1 ' 1  c1.1 1  +  XV"  c'11 
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where  (after  some  algebra) 


T  K'„' 1  Au 


-  KV 

- - — r  [-Xu  •  (t  ©  n) 

I  *:-Xi  I 


ff  +  ff»(t©n)«Au 


+  -r - - - j- Au  •  (n  ©  n)  •  ?t  l(i) 

I  -  x,  | 

xSs) 

VT  K(/'  Au  =  — - - - —  [  I70  •  (n  ©  n)  •  Au  +  lj  •  (n  ©  n)  •  Au 

I  x2  -  X,  | 


I  x: 


ft  •  (n  ©  t  +  t  ©  n)  •  Au  ](J) 


(2.20) 


(2.21) 


2.2.1.  Matrix  representation.  To  facilitate  finite  element  implementation  of  the  above 
derived  tangent  stiffness  operators  and  residuals,  matrix  formulations  are  given  as  follows.  For 
simplicity,  we  will  drop  the  superscript  (i)  and  focus  on  a  typical  single  three-node  element. 


Let  us  start  by  introducing  the  following  vectors 

N®  =  [n  ,  - (1  -  a°)  n  ,  -a0  n]r  (2.22a) 

N,  sc,  =  n  =  [  n  ,  -  ( 1  -  a)n,-an]r  (2.22b) 

Ts  =[t,-(l  -a)t,-a  t]r  (2.22c) 

T  =  7 s  [0  ,  -t ,  t]r  (2.22d) 

N  =  n  =  [ 0  ,  -n  ,  n ]r  (2.22e) 

Au  h  [AUj  ,  Au:  ,  Au2]r  (2.220 


c,  =  t 


(1  -  a°)t- 


gn 


I  x:  -  x,  | 


a°  t  + 


gn 


I  Xf  -  X,  I 


(2.22g) 


where  the  unified  order  of  components  in  all  vectors  has  been:  slave  node  -  master  node  1 
master  node  2.  By  using  these  matrix  notations,  Eq.  (2.20)  and  (2.21)  can  be  rephrased  as 


K 


n 


[NT,r  +  Tj  Nr 


N  Nr] 


K, 


[N"  Nr  +  N  Nsr  - 


r— ^ - — r  (NTrtT  N r ) ] 

I  X;  -  X,  | 


(2.23) 

(2.24) 


In  addition,  from  Eq.  (2.19)  together  with  the  fact  that  -\„  =  u)„  gn  and  =  u.\,  g,  .  the  con¬ 
tact  residual  vector  (due  to  contact  only)  for  a  single  element  is 

Rr  =  -  [<*>„  gn  Cn  -  u),  g,  c.  ]  (2.25) 


The  linearization  of  Eq.  (2.25)  with  respect  to  u  in  conjunction  with  Eq.  (2.19)  then  leads  to 
the  following  perturbed  Lagrangtan  contact  tangent  stitfness  matrix 
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N  N  *  \  N 


I  V  -  V  I 


N  I  ‘  -  T  Nr)l 


Here,  ^  ■  penalty  parameter t  has  been  assumed.  Finally,  the  total  tangent  stiffness 

matrix  and  residual  sector  of  the  bodies  in  contact  take  the  form 


K  =  k*  *  14  K'r 1 

j  -  i 


R  =  -  l  R#  *  ^  *’f?,  c,  *  g,  c,  )lsl  ] 


(2.27) 


(2.28) 


Remark  2.1.  From  Eq.  (2  23).  we  observe  that  K„  is  symmetric.  By  contrast,  from 
1 2.24).  k,  land  hence  k,  in  (2.26))  is  non-symmeiric  as  long  as  S°  #  Nj.  That  is,  only  when 
a  =  a  will  k.  be  symmetric.  This  situation  arises  only  when  the  final  numerical  convergence 
is  achieved  for  the  frictional  non-slip  case.  □ 

Remark  2.2.  If  one  drops  the  nonlinear  terms  (K„  ,  K,)  from  (2  26),  then  the  linearized 
theory  is  recovered.  In  other  words,  only  the  rank-one-update  term  w  (Ns  Sf  +  c,  c,r] 
remains.  □ 

2.3.  Frictional  slip 

For  the  case  of  frictional  slide,  use  of  the  Coulomb’s  law  of  friction  renders 
A;  !  =  u  At  where  u  denotes  the  coefficient  of  friction.  Similar  to  the  development  in  Sec. 
2.2.  characteristic  variational  equations  are: 


ouri(u)  -  A;  ouGn  -  u  A;5UG,  =  0 


SA'  (  -  -  A,  *  G  J  =  0 

it 


(2.29a) 


(2.29b) 


in  which  u  is  a  penalty  parameter.  Note  that  in  Eq.  (2  29a)  the  virtual  work  done  by  the  fric¬ 
tional  force  is  always  negative.  Furthermore,  Eq.  (2.29b)  yields  A„  =  ca  G,  as  the  normal 
penalty  contact  force. 

By  taking  the  variations  at  (u  .  A, ).  Eq.  (2.29a.b)  now  read  (see  (2. 1  Sa.b.c )) 


R/j  *  I\  I  NT  c,"  -  u  V,"  c;.1’)  =  0 


h\l  I  -  —  A.  -  G„)  =  0 


(2.30a) 


(2.30b) 


The  consistent  linearization  of  <2  zda.bi  at  <u  ,  A.,  i  then  yields  the  following  expressions  (see 


i  2.  1  M  i: 


k,  -  n  ik  -  k  i  /i  : c„  -  u  c. 


i  ...  i 


J.W.  Ju  .  R.L.  Taylor  &  L.Y.  Cheng 


8 


Rs  +  M  (.V^cif1  -  a  Wl 

5  -  1 

A  „  +  G„ 

CO 


(2.31) 


Here,  K„  is  the  same  as  Eq.  (2.20)  and  (2.23),  whereas  K,  takes  the  following  matrix  form  (for 
a  single  element) 


K,  = 


a  K 

I  x2  -  x,  | 


[  N?  +  N  N[  - 


gn 

I  X2  -  X,  I 


(NTr+T  Nr)] 


(2.32) 


It  is  noted  that  Eq.  (2.32)  can  be  obtained  by  simply  replace  A,  in  Eq.  (2.24)  by  [-  A„],  as  a 

direct  consequence  of  the  Coulomb's  friction  law. 

Remark  2.3.  If  the  surface  coordinate  a  ~  a°  (i.e.,  approximately  unchanged),  then 
N°  =  Ns  and  K,  is  almost  symmetric.  □ 

Since  \„  =  a >  gn,  the  contact  residual  for  one  element  is 

Rc  =  -  u>  gn  [c„  -  n  c,  ]  (2.33) 


The  linearization  of  (2.33)  at  u  together  with  Eq.  (2.31), (2. 32), (2.23)  then  yield  the  contact 
tangent  stiffness  for  frictional  slip: 


Kc  = 


[Ns  Nf  -  nc,  Nfl- 


gn 


I  x2  -  X,  I 


(NTf  +  T,  Nr  + 


gn 


|  X2  -  X,  | 


-T-M--y  ,  [N?  Nr  +  N  Nf  -  -  g\  (N  Tr  +  T  Nr)] 

I  x2  —  X,  |  |  X2  —  X|  | 


N  Nr] 


(2.34) 


Therefore,  the  total  tangent  stiffness  matrix  and  residual  vector  associated  with  the  contacting 
bodies  are 

5 

K  =  Kg  +  M  K£>  (2.35) 

5*1 

5 

R  =  ~[RB  +  M  wgn{cn  -  mc,P]  (2.36) 

t » i 


Remark  2.4.  It  is  seen  from  Eq.  (2.34)  that  the  contact  tangent  stiffness  for  the  case  of 
frictional  sliding  is  always  non-symmetric  due  to  the  nature  of  the  Coulomb’s  friction  law. 
Even  in  the  event  of  linearized  kinematics  (i.e.,  by  neglecting  nonlinear  terms  K„  and  K,),  this 
is  still  the  case  (see  also  Oden  &  Martins  [1985]).  □ 


3.  Numerical  implementation  and  examples 

In  this  section,  implementation  of  the  proposed  formulation  within  the  context  of  the 
finite  element  method  is  described.  Some  numerical  examples  are  also  presented. 
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3.1.  Finite  element  implementation 

The  master-slave'  slide-line  contact  logic  is  employed,  which  features  both  the  "one- 
pass"  and  “two-pass'  algorithms  (see,  e.g„  Hallquist  [1983]).  For  the  planar  three-node  con¬ 
tact  element  under  consideration,  explicit  vector-component  expressions  (6  d.o.f.)  for  the 
notations  defined  in  Eq.  (2.22a-g)  can  be  obtained.  For  example. 


-s 

c 

II 

r 

U 

III 

(1  -a)s 
-(1  -  a)c 

a  s 

-a  c 

-(1  -a°)c  + 


-(1  -a°)s 


- a  c 


-a  s  + 


I  *2-X,  | 

gn 

I  *2-X|  I 

gn 


I  Xj-X,  | 

gn 

1*2-*.  | 


where  s  ,  c  denote  sin0  ,  cost)  ,  respectively  (see  Fig.  1). 


(3.1) 


In  the  spirit  of  operator  splitting  methodology  for  the  Coulomb’s  law  of  friction,  each 
load  (time)  step  is  decomposed  into  two  parts  :  (i)  By  assuming  a  sticking  condition,  a  'stick 
trial '  step  is  first  performed  (similar  to  the  "elastic  trial"  step  in  classical  elastoplasticity).  If 
the  trial  is  successful,  the  contacting  bodies  are  considered  to  be  in  a  state  of  frictional  stick¬ 
ing.  Otherwise,  (ii)  a  'slip  correction'  step  is  performed  (similar  to  the  "plastic  return  map¬ 
ping"  step  in  elastoplasticity)  and  the  bodies  in  contact  are  viewed  to  be  in  a  state  of  frictional 
sliding.  This  operator  split  treatment  separates  the  no-slip  and  slip  conditions  and  renders  the 
transition  from  stick  to  slip  (or  vice  versa)  exactly  the  same  way  as  the  corresponding  case  in 
classical  elastoplasticity.  The  analogy  between  the  Coulomb's  law  of  friction  (for  stick-slip 
contact  problems)  and  yield  criterion  (for  elastoplasticity)  is  noted 


In  all  numerical  examples  that  follow,  standard  Newton  s  method  is  used  for  solution 
procedure.  It  is  emphasized  that  line  search  plays  no  role  in  numerical  simulations  presented 
in  this  section. 


3.2.  Example  1:  frictional  stick 

This  section  is  concerned  with  a  (rigid  or  deformable)  punch  into  an  elastic  foundation 
under  the  circumstance  of  frictional  stick.  See  Figure  2  for  (plane  strain)  finite  element  mesh 
and  dimensions. 

Case  1.  Rigid  punch.  The  material  properties  employed  in  the  computation  are: 
Epunth  =  10*  (assumed  rigid),  rpun<h  =  0  .  E,„und  =  10'  ,  i>.ouna  =  0.3  ,  and  w  «  10"  (penaltv 
value).  The  one-pass  algorithm  is  used  in  this  example.  The  finite  element  solutions  converge 
quadratically  within  4  iterations,  see  Table  I  for  numerical  performance.  The  deformed  mesh 
is  displayed  in  Figure  3.  in  which  the  deformation  is  enlarged  1000  times  the  real  scale  in 
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order  to  fully  see  the  details. 

Table  1.  Residual  &  energy  norms  for  iterates 


Iteration 

1 

2 

3 

4 

Residual 

.245e+4 

.299e+3 

.7 1 9e—  1 

.173e-6 

Energy 

.552e+2 

.151e+0 

,692e-10 

.727e-19 

Case  2.  Elastic  punch.  The  punching  block  is  now  a  deformable  body.  The  material 
properties  involved  in  the  computation  are:  Epunch  =  104  ,  upunch  =  0.3  ,  Efound  =  103  , 
i 'found  =  0.3  ,  and  u>  =  107  (penalty  value).  The  one-pass  algorithm  is  adopted  in  this  example. 
Once  again,  the  finite  element  solutions  converge  quadratically  within  3  iterations;  see  Table  2 
for  numerical  performance.  The  deformed  mesh  is  displayed  in  Figure  4  (to  scale). 

Table  2.  Residual  &  energy  norms  for  iterates 


Iteration 

1  2  3 

Residual 

,245e+4  .528e+2  ,155e-4 

Energy 

.682e+4  ,2l7e-4  . 1 16e-12 

3.3.  Example  2:  frictional  slide 

Attention  is  now  focused  on  the  event  of  frictional  stick-slip  motion.  The  transition 
from  stick  to  slip  (or  vice  versa)  is  accounted  for  in  this  example.  We  once  more  consider  an 
elastic  punch  on  top  of  an  elastic  foundation  made  of  same  materials.  The  finite  element 
mesh  and  boundary  conditions  are  the  same  as  Sec.  3.2  (see  Fig.  2).  Moreover,  the  material 
properties  used  in  the  simulation  are:  E  =  104  ,  v  =  0  ,  and  u>  =  105  ,  jx  =  0.1  (coefficient  of 
friction). 

The  punch  is  first  vertically  loaded  into  the  elastic  foundation,  then  move  horizontally 
to  the  right  by  displacement  controlled  loading  condition  (vertical  loads  still  remain).  During 
the  initial  vertical  loading,  three  bottom  nodes  of  the  punch  are  in  contact  with  the  founda¬ 
tion.  In  particular,  the  two  (outer)  edge  nodes  of  the  punch  undergo  frictional  slip  while  the 
central  node  experiences  frictional  suck.  The  solutions  converge  in  7  iterations  with  a  residual 
norm  less  than  10'5.  Within  the  proposed  formulation  and  implementation,  tangential 
motion  across  element  boundaries  does  not  impose  numerical  difficulties. 


Finite  element  formulation  of  frictional  contact  problems 


11 


Before  the  first  contact  node  (the  rightmost  contact  node)  of  the  punch  reaches  the  right 
edge  of  the  foundation,  a  typical  iteration  count  for  numerical  convergence  is  6  or  7.  After  the 
first  punch  element  begins  to  overhang,  the  contact  area  is  not  constant  and  the  number  of 
contacting  nodes  changes.  At  5%  overhang  (of  the  first  punch  element),  the  convergence  takes 
7  iterations.  See  Fig.  5  for  deformed  configuration  (to  scale).  At  50%  overhang,  8  iterations 
are  taken  before  convergence  is  achieved.  At  80%  overhang,  1 1  iterations  are  recorded.  At 
90%  overhang,  13  iterations  are  required.  Finally,  at  98%  and  100%  overhang,  15  iterations 
are  observed.  See  Figures  6,  7  for  deformed  meshes  (to  scale).  After  the  first  punch  element 
completely  overhangs,  the  solutions  diverge  which  corresponds  to  the  physical  drop-off  pro¬ 
cess  of  the  punch.  For  this  simulation,  it  is  crucial  to  use  the  two-pass  algorithm  for  a  solu¬ 
tion  to  converge.  The  one-pass  algorithm  works  only  before  the  punch  overhangs.  This  exam¬ 
ple  provides  a  severe  test  for  finite  element  formulation  of  frictional  contact  problems. 

To  assess  the  significance  of  the  proposed  consistent  tangent  stiffness,  we  repeat  the 
above  numerical  experiment  by  using  the  linearized  tangent  (i.e.,  employing  only  the  rank- 
one-update  terms).  Before  the  first  punch  element  overhangs,  numerical  convergence  typically 
takes  8  or  9  iterations.  At  5%  overhang,  the  convergence  takes  9  iterations.  At  50%  overhang, 
1 1  iterations  are  observed  before  convergence  is  achieved.  At  80%  overhang,  19  iterations  are 
recorded.  At  90%  overhang,  25  iterations  are  required.  Finally,  at  100%  overhang,  33  itera¬ 
tions  are  observed.  The  significance  of  the  proposed  formulation  versus  linearized  theory  is 
clearly  demonstrated. 

4.  Conclusion 

On  the  basis  of  an  operator  split,  the  proposed  formulation  accommodates  the  frictional 
stick-slip  motion  in  a  variational  framework.  By  a  consistent  linearization  procedure,  explicit 
expressions  for  the  consistent  contact  tangent  stiffness  and  residual  have  been  obtained.  The 
analogy  between  the  proposed  treatment  for  the  stick-slip  motion  and  the  corresponding  treat¬ 
ment  in  classical  elasto-plasticity  is  noteworthy.  In  addition,  for  infinitesimal  deformations  (as 
a  special  case),  the  proposed  formulation  reduces  to  the  linearized  theory  involving  only 
rank-one-update  terms  in  the  tangent  matrices. 

To  illustrate  the  numerical  performance  of  the  proposed  formulation,  some  numerical 
examples  have  been  presented  in  Sec.  3.  The  significant  role  of  the  proposed  tangent  stiffness 
is  fully  demonstrated. 
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Figure  I.  Geometry  and  definition  of  a  typical  three-node  contact  ele- 
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